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ABSTRACT 

A  finite  element  program  based  on  the  plane  stress 
assumption  is  developed  and  applied  to  elasto-plastic  frac¬ 
ture  problems  involving  monotonical ly  increasing  loads.  The 
program  directly  predicts  the  initiation  and  propagation  of 
fracture  in  the  structure.  That  is,  the  concept  of  stress 
intensity  factor  is  not  utilized  in  the  present  approach. 

The  approach  uses  a  piecewise  linear  approximation  of  the 
actual  stress-strain  curve  for  the  material,  and  the  maximum 
strain  criteria  to  predict  both  the  yield  and  fracture.  An 
incremental  loading  technique  is  employed  to  load  the  struc- 
ture,  and  a  "zero  modulus-unload  reload"’  scheme  is  developed 
to  handle  the  response  of  the  structure  at  fracture.  Com¬ 
parisons  with  published  data  on  a  cracked  panel,  and  the  ex¬ 
perimental  data  obtained  during  this  study  on  tensile  and 
cracked  specimens  show  that  the  finite  element  program  de¬ 
veloped  herein  can  accurately  predict  load  and  deflection  at 
fracture,  load-deflection  curves,  fracture  initiation  loca¬ 
tions,  and  stable  or  unstable  crack  propagation.  This  ap¬ 
proach  is  shown  to  be  highly  dependent  on  the  mesh  density 
in  areas  of  high  strain  gradients. 
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IT  WOULD  BE  HARD  TO  IMAGINE  A  TIME  WHEN 
SOMEONE,  SOMEWHERE,  DID  NOT  LOOK  DOWN  ON 
SOME  BROKEN  SOMETHING  AND  WONDER  WHY. 


FRACTURE  PREDICTION  IN  PLANE  ELASTO- 
PLASTIC  PROBLEMS  BY  THE  FINITE 
ELEMENT  METHOD 

CHAPTER  I 

INTRODUCTION 


1.1  Historical  Background 

On  the  10th  of  January  1954,  a  B.O.A.C.  Comet,  the 
first  pressurized  commercial  jet  airliner,  fell  into  the  sea 
near  Rome,  killing  29  passengers  and  its  crew  of  6.  The 
Comet  fleet  was  immediately  grounded  and  carefully  inspected 
and  modified.  Comet  service  resumed  on  the  23rd  of  March, 
and  16  days  later,  a  second  Comet  was  lost.  Subsequent  in¬ 
vestigations  determined  that  the  Comet  had  met  airworthiness 
requirements  effective  at  the  time;  however,  these  require¬ 
ments,  based  on  static  analysis  and  testing,  were  insuffi¬ 
cient  to  predict  the  type  of  cyclic  failure  experienced  by 
the  Comets. 

As  a  result  of  these  accidents,  fatigue  analysis  and 


testing  became  an  integral  part  of  aircraft  design.  Fatigue 
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analysis  frequently  took  the  form  of  a  damage  accumulation 
theory  such  as  Miner's  rule.  With  Miner's  rule,  the  number 
of  cycles  at  each  stress  level  is  divided  by  the  number  of 
cycles  to  failure  at  that  level.  These  fractions  are  then 
summed  with  failure  indicated  by  a  total  accumulation  of 
one.  Fracture  testing  was  accomplished  by  subjecting  an 
airframe  to  "blocks"  of  loading  which  would  simulate  those 
anticipated  in  actual  service.  The  inaccuracies  in  this 
type  of  fatigue  approach  required  considerable  conservatism 
(along  with  the  associated  high  cost)  to  insure  the  integ¬ 
rity  of  the  structure.  Typically,  airframes  were  required 
to  withstand  one  and  a  half  times  the  maximum  static  load 
and  four  times  the  number  of  cyclic  loads  expected  to  be 
encountered  in  service.  Additionally,  most  fatigue  phil¬ 
osophies  dictated  that  any  cracking  was  to  be  considered  a 
failure . 

The  loss  of  a  U.S.  Air  Force  F-lll  in  1969  initiated 
a  rethinking  of  airframe  design  and  analysis  concepts. ^ 
Failure  in  this  aircraft  was  traced  to  a  small  manufactur¬ 
ing  flaw  in  a  wing  pivot  fitting,  not  to  a  design  induced 
fatigue.  In  a  fashion  reminiscent  of  the  Comet  incidents, 
it  became  apparent  that  static  and  fatigue  concepts  alone 
would  not  predict  the  type  of  failure  incurred  by  this  air¬ 
craft.  It  also  became  clear  that  a  more  efficient  approach 
would  be  to  design  a  structure  to  be  crack  tolerant,  and 
that  some  analytic  method  would  be  necessary  to  accomplish 
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this  goal. 


1.2  Fracture  Mechanics 

Fracture  mechanics  theory  which  had  been  successfully 
applied  to  crack  instability  research  was  adapted  to  this 
new  design  role.  A  popular  form  of  the  fracture  mechanics 
"law”  is 

da/dN  =  C  f (K) 

where  da/dN  is  the  crack  growth  rate  (a  being  the  crack 
length,  N  the  number  of  cycles),  C  is  a  material  constant, 
and  K  is  the  stress  intensity  factor  which  relates  the 
stress  conditions  with  the  crack  length.2  When  coupled  with 
a  damage  tolerance  approach  to  design,  an  initial  flaw  size 
(usually  the  minimum  crack  size  which  can  be  repeatedly  de¬ 
tected  by  nondestructive  inspection)  is  assumed.  The  crack 
is  then  grown  according  to  the  appropriate  fracture  mechan¬ 
ics  "law"  until  it  reaches  a  critical  crack  length.  This 
information  is  then  used  to  establish  inspection  intervals 
for  the  structure.  During  the  course  of  the  design,  any 
area  found  to  have  unacceptably  rapid  crack  growth  (i.e., 
inspection  intervals  are  too  close)  must  be  redesigned. 

The  advantage  of  such  a  procedure  is  that  accurate  analysis 
and  inspection  can  safely  extend  the  structure  to  its  full 
useful  life.  However,  due  to  the  randomness  of  possible 
initiation  sites,  each  area  of  the  structure  must  be  anal¬ 
yzed. 
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Unfortunately,  the  application  of  a  fracture  mechanics 
"law"  to  all  areas  of  a  large  structure  is  a  difficult,  if 
not  impossible,  bookkeeping  task.  The  stress  intensity  fac¬ 
tors  must  be  determined  for  each  area  of  the  structure  for 
each  type  of  loading,  a  time  consuming  effort  even  for  simple 
geometries  and  loadings.  Then  these  factors  must  be  combined 
in  the  proper  manner  to  establish  the  crack  growth  rates. 

O 

Additionally,  Boyd*3  pointed  out  that  the  assumptions  asso¬ 
ciated  with  fracture  mechanics  "laws”  are  more  restrictive 
than  is  generally  realized.  The  most  significant  defect  in 
the  application  of  these  fracture  mechanics  approaches  to 
practical  aircraft  structures  lies  in  their  elastic  formu¬ 
lation.  Structural  metals,  however,  exhibit  a  large  degree 
of  plastic  deformation  ahead  of  the  crack  tip  which  signifi¬ 
cantly  affects  their  response.  Some  attempts  have  been  made 
to  provide  correction  factors  for  this  effect,  but  they  only 
further  the  gap  between  the  physical  phenomenon  and  the  an¬ 
alytical  technique. 

1.3  The  Finite  Element  Method 

If  the  finite  element  method,  which  has  enjoyed  enor¬ 
mous  success  in  the  application  to  structural  mechanics 
problems  over  the  last  three  decades,  could  be  employed 
successfully  to  directly  predict  fracture  in  structures, 
the  shortcomings  of  present  methods  could  be  overcome .  The 
method  has  already  been  used  in  fracture  mechanics  to  cal- 
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culate  stress  intensity  factors.  However,  a  direct  proce¬ 
dure  of  predicting  fracture  would  eliminate  the  laborious 
task  of  calculating  stress  intensity  factors  for  all  regions 
of  a  structure  and  using  them  to  predict  fracture.  Additio¬ 
nally,  the  fracture  mechanics  approach  based  on  the  stress 
intensity  factor  does  not  take  into  account  the  plastic  de¬ 
formation  associated  with  the  fracture  phenomena  in  non- 
brittle  materials.  Since  the  finite  element  method  can  also 
be  used  to  analyze  structural  problems  involving  material 
nonlinearities,  it  seems  obvious  to  employ  the  method  to  di¬ 
rectly  predict  fracture  in  elasto-plastic  materials  using  a 
realistic  stress-strain  law. 

In  the  finite  element  method  (FEM),  a  given  structure 
is  divided  into  substructures,  called  finite  elements. 

These  elements  can  be  of  different  shapes  and  sizes  (a  phy¬ 
sical  continuum  can  be  viewed  as  a  collection  of  smaller 
elements) .  A  typical  element  is  isolated  from  the  collec¬ 
tion  and  its  physical  properties,  such  as  the  stiffness  co¬ 
efficients,  are  developed  using  piecewise  approximation  of 
the  variables.  Then  the  discrete  set  of  equations  govern¬ 
ing  the  complete  structure  are  obtained  by  putting  the  ele¬ 
ment  equations  together.  Generally,  the  accuracy  of  the 
predicted  structural  response  improves  with  the  number  of 
elements  (i.e.,  with  the  decrease  of  element  size),  and  the 
order  of  approximation  which  is  used  to  represent  the  solu¬ 
tion.  With  respect  to  fracture  studies,  the  FEM  offers  a 
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unique  opportunity  to  include  the  effects  of  the  material 
nonlinearities.  The  inherent  flexibility  and  ease  of  ap¬ 
plication  of  the  FEM  suggests  that  it  is  a  valuable  design 
tool  for  direct  prediction  of  fracture  as  well  as  the  res¬ 
ponse  of  structures  in  the  presence  of  cracks. 

1.4  Brief  Review  of  Pertinent  Literature 

The  finite  element  method  has  already  been  used  in 
4—1  1 

many  studies  to  investigate  fracture  processes.  Most  of 
these  investigations  have  centered  on  examining  the  local¬ 
ized  effects  of  fracture  and  the  effects  of  cracks  on  struc¬ 
tures,  rather  than  on  predicting  actual  catastrophic  fail¬ 
ure  loads  and  deflections.  Newman10  studied  the  effects  of 
various  parameters  such  as  the  mesh  size,  strain  hardening, 
and  critical  strain  on  finite  element  fracture  prediction; 
however,  no  experimental  results  were  used  for  comparison. 

On  the  other  hand.  Miller  et  al.11  presented  a  finite  ele¬ 
ment  solution  and  experimental  results  for  a  cracked  panel 
under  monotonically  increasing  stress.  Unfortunately,  the 
finite  element  predictions  did  not  show  close  agreement  with 
experimental  data.  These  predictions  also  varied  signifi¬ 
cantly  with  the  method  of  load  redistribution  at  fracture. 
Furthermore,  the  nodal  uncoupling  method  used  by  Newman, 
Miller  and  others  is  restricted  to  fracture  prediction  along 
lines  of  symmetry. 
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I . 5  Objectives  of  the  Present  Stud\ 


The  present  investigation  is  concerned  with  the  de¬ 
velopment  of  a  finite  element  program  to  directly  predict 
fracture  in  non-brittle  materials  under  monotonically  in¬ 
creasing  loads.  The  procedure  involves  the  use  of  a  piece- 
wise  linearized  stress-strain  curve,  with  an  incremental 
loading.  This  study  also  involved  experimental  investiga¬ 
tion  of  fracture  to  determine  the  accuracy  of  the  numerical 
predictions.  Thus,  the  goal  of  this  research  was  to  deter¬ 
mine  if  the  finite  element  method  could  provide  an  accurate 
and  useable  design  tool  for  the  analysis  of  fracture  and 
crack  growth  in  practical  structures.  To  accomplish  this 
goal,  simple,  non-trivial  two  dimensional  (plane  stress) 
structures  under  uniaxial  loading  were  considered.  This 
study  further  defines  the  factors  affecting  accurate  pre¬ 
diction  of  fracture  by  the  finite  element  method,  and  de¬ 
termines  if  the  results  obtained  by  Miller  represent  typi¬ 
cal  errors  to  be  expected  by  such  a  method.  Accuracy  of 
the  method  is  demonstrated  by  comparison  with  experimental 
data . 


CHAPTER  II 


THEORETICAL  CONSIDERATIONS 

II. 1  Governing  Equations 

In  a  continuum,  application  of  loads  results  in 
stresses.  At  any  point  in  the  structure,  there  are  nine 
stress  components;  however,  only  six  of  them,  three  normal 
stresses  (ax>  ay,  az)  and  three  shearing  stresses  (aXy, 
ayz,  azx) ,  are  independent.  The  stresses  induce  strains  in 
the  material.  For  a  three  dimensional  linear  elastic  aniso¬ 
tropic  material,  the  six  strains  (Ej.)  are  related  to  these 
stresses  as  follows: 

Ei  =  Si j  <7j  +  QfT  (2.1) 

where  Sjj's  are  the  compliance  coefficients,  di's  are  the 
coefficients  of  thermal  expansion,  and  T  is  the  temperature 
change.  There  are  36  compliance  coefficients,  but  due  to 
symmetry  of  Sjj  only  21  are  independent.  For  isotropic 
materials  the  number  of  independent  coefficients  is  two. 

Here  it  is  assumed  that  the  material  is  isotropic  and 
the  temperature  changes  are  negligible.  Since  only  thin 
sections  are  to  be  modeled,  a  state  of  plane  stress  (with 
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the  stresses  and  strains  in  these  regions.  In  the  present 
study,  where  aluminum  (2024-T3)  was  used,  it  is  assumed  that 
the  nonlinear  elastic  portion  is  negligibly  small. 

These  two  remaining  regions  are  shown  in  Figure  2.1 
for  a  uniaxial  stress  state.  In  the  initial  linear  region. 


the  material  response  is  elastic  and  structures  whose  loads 
result  in  stresses  in  this  region  will  return  to  their  ori¬ 
ginal  shape  when  the  loads  are  removed.  Structures  loaded 
into  the  plastic  zone,  however,  take  on  a  permanent  set  on 
unloading  (dashed  line  of  Figure  2.1). 

In  order  to  analyze  the  nonlinearity  introduced  by 
the  plastic  response,  the  curve  in  Figure  2.1a  is  divided 
into  a  series  of  linear  portions  as  shown  in  Figure  2.1b, 
with  the  tangent  modulus  and  incremental  Poisson's  ratio 
replacing  the  elastic  constants  previously  mentioned. 

Next,  the  choice  of  failure  criteria  to  determine 
yield  and  fracture  should  be  considered.  The  commonly 
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used  failure  criterion  is  that  of  von  Mises,  which  shows  the 
best  agreement  with  experimental  yield  data  for  metals. 

For  the  case  of  plane  stress,  von  Mises  criterion  takes  the 


form 


ai2  -<*io2  +  o22  =  oj 


(2.6) 


where  a-^  and  02  are  the  principal  stresses,  and  ayp  is  the 
yield  stress.  If  the  left  hand  side  is  less  than  <7  2  yield 

yp 

does  not  occur.  The  surface  described  by  Eq.  (2.6)  is  shown 
in  Figure  2.2.  Problems  arise  in  extending  von  Mises'  cri¬ 
terion  into  the  plastic  range  with  the  incremental  approach 
used  in  this  study  (this  will  be  discussed  in  detail  in 
Chapter  III);  therefore,  the  maximum  strain  criterion  was 
used.  In  this  theory,  yield  occurs  when  the  maximum  strain 
exceeds  the  strain  at  yield.  That  is 


E  -1  =  +  C 

1  -  yp 

C2  =  —  Eyp 
c3  =  i  cyp 


(2.7) 


where  E-^,  £2,  and  e3  are  the  principal  strains.  Conversion 
of  these  equations  to  equivalent  principal  stresses  is  also 
shown  in  Figure  2.2.  As  can  be  seen,  when  0i»02  or 
02»Oi,  both  theories  give  approximately  the  same  results. 
Since  the  stress  fields  in  the  parts  to  be  analyzed  meet 
this  requirement,  the  use  of  the  maximum  strain  criterion  is 
justified  for  this  study. 

The  maximum  strain  criterion  can  also  be  extended 
along  the  stress-strain  curve  to  predict  subsequent  changes 
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in  modulus  and  finally  the  fracture. 

The  principal  strains  used  in  Eq .  2.7  can  be  obtained 
from  the  following  equation: 


1,2 


'x+ey 


a 


£x~£ 


Z>2 


(eXy)‘ 


(2.8) 


A  similar  equation  may  be  used  to  obtain  principal  stresses 
With  the  basic  continuum  equations  developed  here,  a 
suitable  procedure  must  be  employed  to  obtain  a  solution. 


II. 2  The  Finite  Element  Approach 

The  finite  element  method  (FEM)  is  employed  to  solve 
the  elasticity  equations  for  each  increment  of  load.  Only 
a  brief  discussion  of  the  method  will  be  presented  here; 
however,  for  a  more  thorough  presentation,  see  Pef.  12. 

The  basic  element  used  in  this  study  is  the  standard 
constant  strain  triangle  (CST)  shown  in  Figure  2.3.  The 
element  has  a  total  of  six  displacement  degrees  of  freedom, 
one  in  each  direction  at  each  of  three  nodes.  The  displace¬ 
ment  field  is  approximated  by  the  linear  relation  of  the 
form 

u=  £ui4*i  (x,y),  v  =X!vi^i(x>y)  (2-9) 

where  u^,  v^  are  the  nodal  values  of  the  deflections  (at 
node  i),  and  the  M^'s  are  element  interpolation  functions, 
given  by 

>fi  -  5S(ai  +  +  V> 


(2.10) 
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where  A=  area  of  the  triangle, 

ai=  xjyk  "  xkyi’  Pi  =  yj  -  yk>  Yi  =  xk  -  xj .  and 
(Xi.yi)  are  the  coordinates  of  node  i. 

Combining  Eqs.  (2.4)  and  (2.9)  gives: 
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dv  +  Force  Terms 


and  substituting  from  Eq.  2.13, 
1  = 


dv  +  Force  Terms 


(2.16) 


(2.17) 


Applying  variational  methods  to  this  equation  to  minimize 
I  (81=0)  results  in  the  governing  equation  for  this  system. 

Jf|  =  [k] J u  |  (2.18) 

where  |f|  is  the  force  vector,  |^k]  is  the  stiffness  matrix, 
and  juj  is  the  displacement  vector.  The  form  of  £k]  results 
from  the  variation  of  Eq. (2.17) and  is  given  as 

[K]  =  j  i[b]T[e][b]  dv  (2.19) 

which  for  the  CST  element  becomes 

[k]  =  [b]T[E][b]  At  (2.20) 

where  A  is  the  area  and  t  is  the  thickness. 

The  procedure  then  becomes  to  assemble  the  stiffness 
matrix  [k],  element  by  element.  These  are  used  to  assemble 
the  global  stiffness  matrix  for  the  entire  system.  Boundary- 
conditions  must  be  applied  to  the  assembled  system  of  equa¬ 
tions  of  the  form  of  Eq.  (2.18)  before  these  equations  are 
solved  for  juj.  If  strain  or  stress  values  are  desired 
Eq.  (2.12)  and  Eq.  (2.3)  may  then  be  applied. 

This  procedure  is  then  automated,  and  the  analyst  only 
need  describe  the  geometry  in  terms  of  elements  and  nodes 
and  boundary  conditions. 

The  programs  that  incorporate  this  development  are 


described  in  the  next  chapter. 


CHAPTER  III 


FINITE  ELEMENT  FRACTURE  PROGRAMS 

III.l  Introduction 

Three  two-dimensional  plane  stress  finite  element  pro¬ 
grams  are  developed  herein  to  predict  yield  and  fracture 
under  monotonically  increasing  loads.  These  programs  are: 

1.  FRACTURE:  This  finite  element  program  is  devel¬ 
oped  to  analyze  point  loaded  tensile  and  notched  specimens. 
Engineering  stress-strain  relations  are  used;  however,  the 
model  geometry  is  not  updated  during  each  load  increment . 

2.  PANEL1:  This  program  is  a  modification  of  FRAC¬ 

TURE  and  is  used  to  analyze  uniformly  loaded  panel  specimens 

3.  PANEL 2 :  This  is  a  modification  of  PANELl  which 

uses  incremental  geometry  changes  and  true  stress-strain 
relationships . 

Input  was  obtained  from  a  mesh  generation  program, 
and  all  input  data  was  plotted  as  a  check  for  errors.  All 
of  the  programs  were  run  on  the  University  of  Oklahoma's 

Merrick  Center  IBM  370/158  computer. 
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The  programs  contained  in  this  work  are  not  optimized 
or  even  necessarily  efficient  from  the  programming  point  of 
view. 

III. 2  Formulation 

The  basis  of  the  formulation  of  the  programs  developed 
herein  is  that  each  element  of  the  finite  element  mesh  has 
its  own  material  properties  (modulus  and  Poisson's  ratio) 
based  on  its  state  of  strain.  These  "local"  properties 
should  approximate  those  of  the  actual  structure.  Further¬ 
more,  these  properties  will  be  those  of  a  uniaxial  tensile 
test  specimen  of  the  same  material  under  the  same  state  of 
strain;  that  is,  when  an  element  has  a  principal  strain 
equal  to  the  uniaxial  yield  strain  of  the  material,  the  ele¬ 
ment  yields  (changes  tangent  modulus  and  Poisson's  ratio). 

In  a  similar  manner,  an  element  fractures  (changes  modulus 
to  zero)  when  its  maximum  principal  strain  is  equal  to  the 
strain  at  which  the  tensile  specimen  fractures. 

To  apply  these  concepts  to  an  operational  program,  it 
is  necessary  to  have  the  entire  stress-strain  curve  from 
the  elastic  region  all  the  way  to  fracture  (while  stress- 
strain  curves  are  readily  available,  strain  at  fracture  is 
not)  .  The  stress-strain  curve  is  then  divided  into  a  series 
of  linear  segments  as  shown  in  Figure  3.1.  From  this  lin¬ 
earized  curve,  values  for  modulus  are  obtained  as 
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gj  ~  o i-1 
ci  "  Ei-1 


(3.1) 


where  a±  and  are  the  engineering  stress  and  strain,  res¬ 
pectively  . 

The  values  of  Poisson's  ratio  (  jj  )  in  the  plastic 

range  were  calculated  from  the  equation  given  by  Bert,  Mills, 
1 3 

and  Hyler  as  proposed  by  Nadai:^ 

M  =  h  -  (i  -kfe)(Es/Ee)  (3.2) 

where  Me  is  the  elastic  Poisson's  ratio,  Eg  is  the  secant 
modulus,  and  Ee  is  the  elastic  modulus.  For  any  section, 
the  secant  modulus  is  taken  as  the  average  stress  in  the 
interval  divided  by  the  average  strain  in  that  interval; 
that  is. 


Therefore, 


ES.  =  qi  +  qi-l 
€i  +  C  i-1 

the  incremental  Poisson's  ratio  (MO 
u  •  =  i  _  ~Me)  (Oj  +  Oj-l) 

1  Ee  ( Ci  +  Gj^_j) 


(3.3) 


becomes 


(3.4) 


These  material  properties  along  with  the  terminal 
strains  (C^)  for  each  interval  are  stored  in  the  program 
and  referenced  by  an  element  material  pointer.  The  proper¬ 
ties  for  a  fractured  element  are  also  stored,  with  the  tan¬ 
gent  modulus  set  to  zero  and  Poisson's  ratio  equal  to  0.5. 
As  an  example  Table  3.1  gives  the  tabulated  steel  and  alum¬ 
inum  properties  used  in  the  program  FRACTURE. 

With  the  material  properties  tabulated  and  referenced 
by  strain  level,  the  material  nonlinearity  of  the  problem 
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TABLE  3.1 


EXAMPLE  OF  TABULATED  MATERIAL  PROPERTIES 


Pointer 

Region 

Ei 

Mi 

Ci 

1,1* 

Elastic 

29000000. 

0.318 

0.00200 

1,2 

Elastic 

10500000. 

0.313 

0.00472 

2,2 

1st  yield 

1183000. 

0.375 

0.00753 

3,2 

2nd  yield 

462000. 

0.443 

0.01814 

4,2 

3rd  yield 

271000. 

0.467 

0.03286 

5,2 

4th  yield 

152000. 

0.484 

0.07817 

6,2 

5th  yield 

65000. 

0.491 

0.14291 

7,2 

Fracture 

0. 

0.500 

999. 

♦First  #=Region,  Second  #=Material  (l=Steel,  2=Aluminum) 


is  approximated  by  considering  the  structure  to  be  analyzed 
as  a  composite  of  a  finite  number  (n)  of  elements  with  ap¬ 
propriate  properties.  At  the  start  of  the  analysis,  the 
entire  composite  is  assumed  to  have  the  same  properties 
(those  of  the  elastic  portion  of  the  stress-strain  curve). 

As  the  load  increases,  one  element  (say  the  k-th  element) 
will  reach  a  total  principal  strain  value  equal  to  the  yield 
strain.  This  kth  element's  properties  (modulus  and  Pois¬ 
son's  ratio)  are  modified;  therefore,  the  composition  of 
the  structure  is  n-1  elastic  elements  and  one  element  with 
a  reduced  (plastic)  modulus.  Since  the  response  of  each 
element  is  again  linear,  the  usual  elastic  finite  element 


22 


analysis  can  be  used  until  another  element  yields  or  the  kth 
element  reaches  the  strain  at  which  mext  modulus  change 
occurs.  The  process  is  repeated  until  at  least  one  element 
fractures.  The  total  load  to  this  point  is  the  sum  of  the 
incrementally  applied  loads,  and  the  total  deflections  at 
any  node  are  the  sum  of  the  incremental  deflections. 

At  fracture,  the  procedure  described  above  must  be 
modified,  since  the  loads  carried  by  the  fractured  element 
to  this  point  must  now  be  carried  by  the  remaining  structure 
In  Ref.  10  anc!  11  the  crack  is  advanced  by  removing  the  con¬ 
straint  at  the  fractured  node  and  redistributing  the  force 
on  that  node  to  the  remaining  nodes  along  the  appropriate 
line  of  symmetry;  however,  as  Miller11  points  out,  there  is 
no  obvious  rationale  which  appears  to  govern  the  redistri¬ 
bution,  while  the  effects  of  the  redistribution  method  are 
significant.  In  the  present  program,  a  new  and  completely 
general  procedure  is  applied.  When  an  element  reaches  the 
strain  for  transition  to  fracture,  the  structure  is  unloaded 
following  the  elastic  response  of  the  unfractured  specimen 
while  retaining  each  element's  progress  along  the  stress- 
strain  curve.  It  is  then  reloaded  with  the  fracture  sur¬ 
face  extended.  Any  effects  due  to  compression  during  this 
unloading  are  ignored  since  the  actual  structure  never  under¬ 
goes  this  unload  reload  cycle.  If  the  main  diagonal  stiff¬ 
ness  coefficient  corresponding  to  a  node  is  reduced  to  zero 
(the  node  is  unconnected)  as  a  result  of  a  modulus  change 


at  fracture,  the  node  is  constrained. 

If  another  element  fractures  before  the  maximum  load 
is  reachieved,  then  the  crack  growth  is  unstable  but  may  be¬ 
come  stable  again  if  the  load  subsequently  increases  over 
the  previous  maximum  fracture  load.  Figure  3.2  illustrates 
this  procedure  for  a  stress-strain  curve  with  3  linear  plas¬ 
tic  regions.  A  sample  mesh  at  a  crack  tip  is  shown  in 
Figure  3.2  along  with  four  sample  elements  numbered.  The 
bottom  four  curves  are  plots  of  typical  stress-strain  res¬ 
ponses  of  the  four  elements  as  load  increases.  Numbers 
along  the  curves  indicate  the  load  level  at  that  point.  At 
level  1,  element  1  at  the  crack  tip  enters  the  first  yield 
region.  There  is  less  stress  (strain)  concentration  at  the 
other  three  elements;  therefore,  they  advance  only  partially 
along  the  elastic  portion  of  the  curve.  At  load  level  2, 
element  1  changes  modulus  again,  even  before  any  of  the 
other  three  elements  have  yielded  for  the  first  time.  The 
load  continues  to  increase  to  level  3,  at  which  point  the 
principal  strain  of  element  2  indicates  it  has  reached  the 
transition  strain  for  first  yield.  This  continues  through 
levels  4,  5,  and  6.  At  level  7,  element  1  has  reached  the 
fracture  strain.  The  entire  model  is  then  unloaded  (arti¬ 
ficially).  The  modulus  for  element  1  is  set  to  zero  and 
reloading  begins  with  all  unfractured  elements  having  an 
elastic  modulus.  Stress  remains  at  the  unloaded  value  for 
element  1,  since  all  incremental  stresses  are  zero.  Strains 
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for  element  1,  however,  continue  to  increase.  At  level  8 
element  4  changes  modulus  for  the  first  time.  At  level  9, 
the  stress  and  strain,  but  not  the  load,  at  element  2  equals 
the  equivalent  values  at  level  7.  Here  the  element's  prop¬ 
erties,  just  at  the  elastic  values,  return  to  those  of  yield 
zone  2.  Elements  3  and  4  retain  the  elastic  modulus  and 
Poisson's  ratio.  At  level  10,  element  2  enters  the  3rd 
yield  zone  and  at  11,  element  4  enters  the  2nd  yield  zone. 

At  level  12,  element  2  fractures  and  the  structure  is  again 
artificially  unloaded.  Note  that  element  3  properties  have 
remained  elastic  since  element  1  fractured.  At  13,  reload¬ 
ing  has  begun  but  at  a  slower  rate  for  element  3  due  to  the 
low  strain  behind  the  crack  tip.  If  the  load  at  level  12 
is  larger  than  the  load  at  7,  then  fracture  at  7  is  stable. 
On  the  other  hand,  if  load  7  exceeds  level  12,  then  fracture 
is  unstable.  This  process  continues  until  the  stiffness 
matrix  is  no  longer  invertable. 

The  "zero  modulus-unload  reload"  method  just  described 
has  the  following  advantages  over  the  nodal  release-load 
redistribution  approach  of  Ref.  10  and  11: 

1.  The  method  of  redistribution  is  not  arbitrary, 
but  based  on  cracked  specimen  geometry. 

2.  Failures  can  occur  anywhere  in  the  model.  With 
sufficiently  small  elements,  the  zero  modulus  elements  act 
as  the  crack.  Not  only  is  the  nodal  release  method  incap¬ 
able  of  predicting  general  crack  growth,  it  is  specifically 
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restricted  to  the  study  of  fracture  along  axes  of  symmetry. 

3.  The  "zero  modulus-unload  reload"  method  described 
herein  is  element  oriented,  whereas  the  load  redistribution 
method  is  node  oriented.  Since  only  deflections  are  speci¬ 
fied  at  the  nodes  and  stresses  and  strains  are  specified 
over  the  elements,  stress  or  strain  data  must  be  arbitrarily 
distributed  to  the  nodes  in  the  load  redistribution  method 
so  that  the  failure  criteria  may  be  applied. 

The  current  method  does  require  the  extra  time  used 
to  unload  and  reload  the  structure. 

This  study  uses  maximum  strain  criteria  both  for 
yielding  and  fracture.  Newman1-0  and  Miller  et  al  A1  both  use 
maximum  strain  for  fracture,  but  use  von  Mises  criteria  for 
yielding.  A  problem  with  the  stress  formulation  of  von  Mises 
criteria  is  that  not  all  practical  materials  (  for  example 
mild  steel)  have  unique  strains  for  a  given  stress.  Thus 
if  the  material  stress-strain  curve  is  as  shown  in  Figure  3.3, 
and  one  linearized  section  is  taken  from  a  to  b,  then  the 
modulus  for  section  ab  is  zero.  The  incremental  stress  for 
an  element  with  properties  in  this  portion  of  the  curve  will 
always  be  zero.  Therefore,  under  the  stress  formulation  of 
von  Mises  criteria,  the  stress  will  never  advance  past  point 
b.  This  problem  could  be  circumvented  by  reformulating 
von  Mises  criteria  in  terms  of  strain.  Since  the  specimens 
for  this  study  were  uniaxially  loaded  and  the  minimum  stress 
was  low  compared  to  the  maximum  stress,  maximum  strain 
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criteria  and  von  Mises  criteria  are  very  similar,  and  there¬ 
fore,  this  reformulation  was  not  attempted.  It  should  be 
remembered,  however,  that  for  geometries  and/or  loadings 
which  result  in  approximately  equal  principal  stresses,  the 
principal  strain  criteria  used  in  this  program  will  intro¬ 
duce  significant  error. 

As  described  in  Chapter  II,  the  programs  developed 
here  utilize  a  standard  constant  strain  (linear  deflection) 
triangular  element,  CST,  under  plane  stress  conditions. 

This  element  enables  the  use  of  a  large  number  of  elements 
in  a  given  area  with  minimum  storage  requirements  and  mesh 
refinement  is  easier  to  accomplish.  The  predicted  rate  of 
crack  growth  must  be  independent  of  the  element  size  (crack 
growth  of  0.001  inches  cannot  be  predicted  using  elements 
with  sides  0.1  inches  long);  therefore,  many  small  elements 
are  needed  along  anticipated  crack  paths.  If  higher  order 
elements  (e.g.,  the  linear  strain  triangle)  are  used  in  such 
a  dense  mesh,  storage  requirements  become  excessive.  There¬ 
fore,  the  simplest  two  dimensional  element  (the  CST)  is  used 

In  the  experimental  procedure,  the  tensile  specimens 
are  loaded  through  a  steel  pin  in  the  center  of  the  head  of 
the  specimens.  To  realistically  simulate  this  composite 
structure,  the  finite  element  analysis  includes  the  pin  as 
part  of  the  system.  One  half  of  the  pin  is  divided  into 
six  elements  in  the  one  quadrant  models,  and  a  full  pin  is 
divided  into  twelve  elements  in  the  half  specimen  models. 


A  unit  load  is  applied  to  the  center  of  the  pin  to  distrib¬ 
ute  the  load.  In  programs  PANEL  1  and  PANEL  2,  the  load  is 
uniformly  distributed  to  the  nodes  along  the  edge  of  the 
panel . 


III. 3  Description  of  the  Computer  Programs 

Each  of  the  three  programs  consists  of  a  main  program 
and  six  major  subroutines  as  shown  in  Figure  3.4.  Subrou¬ 
tines  FAIL  and  CHANGE  are  the  only  routines  not  common  to  a 
standard  elastic  finite  element  analysis.  Appendix  I  con¬ 
tains  a  listing  of  FRACTURE  with  significant  differences 
between  FRACTURE,  PANEL 1 ,  and  PANEL2  discussed  in  this 
section . 

The  main  program  first  calls  IREAD,  which  reads  in  the 
geometric  description  (nodal  locations  and  element  conduc¬ 
tivity)  from  the  mesh  generator  as  well  as  the  nodal  con¬ 
straints.  IREAD  also  prints  out  the  data  as  a  check  on 
proper  input . 

The  subroutine  PROP  sets  up  the  material  property  mat¬ 
rix.  This  matrix  contains  the  modulus,  Poisson's  ratio,  and 
strain  for  next  transition  indexed  by  a  pointer  and  type  of 
material .  Recall  that  programs  FRACTURE  and  PANEL  use  mater¬ 
ial  properties  that  are  based  on  the  engineering  stress- 
strain  relation  while  program  PANEL2  uses  the  true  stress- 
strain  relation.  The  subroutines  IREAD  and  PROP  are  called 
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only  once  in  the  program  and  no  further  reference  is  made 
to  them  later  in  the  program. 

The  subroutine  ASSEM  assembles  the  element  stiffness 
properties  to  obtain  the  global  stiffness  matrix,  GSTIF. 

For  each  element,  ASSEM  calls  STIFF  which  calculates  the 
element  stiffness  matrix.  The  global  stiffness  is  stored 
in  a  banded  form  in  the  interest  of  storage  and  computational 
efficiency.  The  appropriate  boundary  conditions  on  the 
nodal  deflections  are  then  applied  by  calling  the  subroutine 
BNDRY .  Finally,  a  check  is  made  to  insure  that  no  main  di¬ 
agonal  terms  are  zero.  This  occurs  if  the  stiffness  con¬ 
tribution  of  each  element  touching  the  node  is  zero;  if 
any  diagonal  term  is  zero,  the  node  is  condensed  out  by 
BNDRY . 


Subroutine  S0LVE1  solves  the  banded  system  of  equa¬ 


tions  , 


-  M  l,»| 


uj.  Here  |f 


(3.4) 
is  the  nodal 


for  the  unknown  displacements 

force  vector  and  is  the  global  stiffness  matrix.  The 
program  was  originally  developed  using  a  Gauss-Seidel  iter¬ 
ative  solution  technique  since  and  juj  change  very  little 

from  load  increment  to  increment.  This  eliminates  the  need 
to  reassemble  the  stiffness  matrix  for  each  load  increment. 
For  a  small  mesh,  the  iterative  method  converged  rapidly 
for  the  first  few  iterations,  but  the  time  required  for  accu¬ 
rate  solutions  greatly  exceeded  that  required  for  the  Gauss 
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elimination  method.  As  a  result,  the  Gauss  elimination 
method  was  used  exclusively  in  this  study;  however,  for  the 
very  large  meshes  used  in  the  following  chapters,  this  iter¬ 
ative  scheme  may  deserve  more  attention.  When  the  specimen 
fails,  the  stiffness  matrix  is  no  longer  invertable  and  the 
computation  is  terminated.  The  load  matrix  for  FRACTURE 
consists  of  a  unit  load  applied  to  the  center  node  of  the 
pin  (used  to  distribute  the  load  to  the  specimen)  in  the  lon¬ 
gitudinal  direction.  For  the  PANEL  programs,  where  the 
specimen  is  uniformly  loaded,  the  unit  load  is  divided  be¬ 
tween  the  top  five  nodes  of  the  specimen  in  ratios  of  1/8, 
1/4,  1/4,  1/4,  and  1/8  starting  from  the  edge  node.  Double 
precision  is  used  in  this  subroutine  and  throughout  the  pro¬ 
gram  to  reduce  roundoff  errors. 

Subroutine  FAIL  is  the  first  non-standard  subroutine 
of  FRACTURE;  i.e.,  it  cannot  be  found  with  the  usual  finite 
element  analysis  programs.  First,  FAIL  calculates  the  unit 
strains  in  the  x  and  y  directions  along  with  the  shear 
strains  for  each  element.  These  are  a  function  of  the  de¬ 
flections  calculated  in  S0LVE1.  The  first  time  through 
FAIL,  all  elements  are  in  the  elastic  range  and  the  initial 
strains  are  zero.  Therefore,  the  principal  strains  for  a 
unit  load  are  calculated  and  the  load  at  first  yield  is 
taken  to  be  the  tabulated  strain  at  first  yield  divided  by 
the  maximum  principal  strain  for  a  unit  load.  Next,  the 
unit  strains  are  multiplied  by  the  calculated  load  to  obtain 


the  total  strains  for  each  element  at  first  yield. 

For  subsequent  calls  to  FAIL,  the  total  strains  for 
each  element  are  not  zero  and  their  directions  are  not  the 
same  as  those  of  the  incremental  strains.  This  makes  di¬ 
rect  calculation  of  the  next  yield  or  fracture  load  impos¬ 
sible;  therefore,  an  incremental  scheme  must  be  used  to 
predict  the  next  load  increment.  This  is  accomplished  by 
storing  the  next  strain  for  modulus  change  for  each  element 
in  an  array  called  ECH ,  and  calculating  the  incremental 
principal  strain  (DEPR)  and  total  principal  strain  (EPRI). 

If  the  incremental  principal  strains  are  in  the  same  direc¬ 
tion  as  the  total  principal  strains,  the  incremental  load 
(LI)  for  the  next  failure  is  given  by 

LI  =  (  ECH  -  EPRI  )/  DEPR  (3.5) 

Since  EPRI  and  DEPR  are  not  in  general  in  the  same  direction, 
Eq.  (3.5)  is  only  an  approximation.  This  calculation  is 
made  for  each  element  and  the  smallest  load  increment  is  then 
used  as  the  trial  load  increment  which  will  cause  the  next 
element  to  change  modulus.  In  order  to  reduce  the  compu¬ 
tational  time,  this  increment  may  be  increased  to  cause 
more  elements  to  fail  for  each  solution  of  Eq .  (3.4).  With 
the  incremental  load  now  calculated,  the  total  strains  are 
set  equal  to  the  previous  total  strains  plus  the  incremental 
strains  times  the  incremental  load.  The  total  principal 
strains  for  each  element  are  then  calculated  and  compared 
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with  the  tabulated  strains  for  next  modulus  change.  If  no 
elements  exceed  the  next  change  strains,  the  incremental 
strains  are  again  multiplied  by  the  incremental  load  and 
added  to  the  total  strains.  The  principal  strains  are  again 
calculated  and  are  compared  to  the  tabulated  values.  The 
cycle  continues  until  at  least  one  value  of  the  tabulated 
maximum  allowable  strain  for  the  interval  is  exceeded.  All 
such  elements  are  printed  along  with  the  new  total  load  and 
the  corresponding  property  intervals.  If  only  yielding  has 
occurred,  control  is  returned  to  the  main  program.  If  frac¬ 
ture  has  occurred,  the  total  strains  are  reduced  by  the  elas¬ 
tic  response  due  to  a  unit  load  times  the  total  load  at  frac¬ 
ture,  and  then  control  is  returned  to  the  main  program. 

Subroutine  CHANG  was  originally  conceived  to  update 
the  global  stiffness  matrix  for  failed  elements  in  conjunc¬ 
tion  with  the  SOLVE  (Gauss-Seidel)  iterative  routine.  Since 
the  entire  stiffness  matrix  must  now  be  regenerated  (because 
it  is  changed  during  the  Gauss  elimination  in  S0LVE1)  for 
each  pass,  the  function  of  CHANG  has  been  reduced  to  up¬ 
dating  the  material  pointer  for  failed  elements.  Element 
thickness  was  also  updated  in  CHANG  for  the  PANEL2  program. 

The  program  continues  to  cycle  from  ASSEM  to  S0LVE1 


to  FAIL  to  CHANG  until  the  stiffness  matrix  can  no  longer 
be  inverted  or  the  program  exceeds  the  estimated  time  limit. 


CHAPTER  IV 


CRACKED  PANEL  ANALYSIS 

IV. 1  Introduction 

After  the  present  study  was  undertaken,  and  the  FRAC¬ 
TURE  program  was  completed,  the  article  by  Miller11  on  frac¬ 
ture  prediction  became  available.  Miller's  work  raised  two 
questions  pertinent  to  the  present  study.  First,  can  the 
present  program  improve  on  Miller's  prediction  of  fracture 
load,  and  second,  can  the  program  accurately  predict  the  ex¬ 
perimental  results  presented  by  Miller? 

To  answer  these  questions  and  to  determine  the  parame¬ 
ters  that  effect  FEM  fracture  prediction,  a  modification  of 
FRACTURE,  called  PANEL1,  was  made.  The  essential  differences 
between  FRACTURE  and  PANEL1,  the  method  of  load  application 
and  the  size  of  the  incremental  load  steps,  are  minor.  A 
third  program,  PANEL2 ,  used  true  stress-strain  relations 
and  updated  geometric  coordinates  and  element  thicknesses 
for  each  load  increment. 
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IV. 2  Finite  Element  Analysis 

Due  to  the  geometric  and  loading  symmetry  of  the  panel, 
only  one  quadrant  of  the  panel  is  used  in  the  finite  element 
analysis.  Figure  4.1a  shows  the  dimensions  of  the  0.10", 
2024-T3  aluminum  test  panel.  A  nonuniform  finite  element 
mesh  of  the  panel  is  shown  in  Figure  4.1b.  To  predict  the 
response  at  the  crack  tip  more  accurately,  further  refine¬ 
ment  was  made  there,  as  shown  in  Figure  4.1c.  The  crack- 
tip  portion  is  blown  up  in  Figure  4. Id  to  show  the  mesh  de¬ 
tails.  The  stress-strain  curve  data  for  the  material  was 
obtained  from  published  data1^.  Some  error  is  introduced  by 
this  selection  since  the  exact  material  properties  are  not 
known.  The  stress-strain  curve  was  then  divided  into  one 
elastic  and  three  linearly  plastic  regions.  The  panel  was 
first  analyzed  using  a  coarse  mesh.  Subsequent  refinements 
of  the  mesh  were  made  at  the  crack  tip.  The  loads  predicted 
by  each  of  these  meshes  are  shown  in  Figure  4.2  as  a  function 
of  the  minimum  element  area  at  the  crack  tip.  Entry  into 
each  plasticity  region  of  the  stress-strain  curve  is  shown 
by  the  lower  three  curves  with  initial  fracture  and  final 
fracture  shown  in  the  top  curves.  The  horizontal  line  repre¬ 
sents  Miller ' s^-  experimental  results.  The  elements  on  the 
right  side  of  the  figure  are  too  large  to  predict  stable 
fracture;  therefore,  only  initial  fracture  is  shown  in  this 
area.  As  can  be  seen  in  the  figure,  as  element  size  becomes 
smaller  at  the  tip,  the  load  at  entry  into  each  of  the 
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plasticity  regions  decreases,  and  the  predicted  final  frac¬ 
ture  curve  converges  to  the  experimental  value.  Note  that 
there  is  no  significant  difference  in  the  present  solution 
and  that  of  Miller's  for  equivalent  mesh  sizes.  The  large 
differences  between  Miller's  experimental  and  numerical  re¬ 
sults  appears  to  be  due  solely  to  the  refinement  of  the  mesh. 
While  the  methods  do  appear  to  have  similar  accuracies,  it 
should  be  remembered  from  earlier  discussions  that  Miller's 
method  only  applies  to  failure  along  lines  of  symmetry,  and 
load  redistribution  procedures  are  arbitrary. 

For  stable  fracture  prediction,  the  mesh  not  only 
needs  to  be  refined  at  the  crack  tip,  but  also  along  the 
projected  crack  path.  The  element  meshes  for  the  data 
shown  in  Figure  4.2  are  basically  the  same  except  for  re¬ 
finements  at  the  crack  tip.  The  mesh  shown  in  Figures 
4.1c  and  4. Id,  however,  is  refined  along  the  entire  path 
of  anticipated  stable  fracture.  While  the  element  size  at 
the  crack  tip  is  larger  for  this  fine  mesh  than  those  at 
the  extreme  left  of  Figure  4.2,  the  predictions  are  more 
accurate  as  can  be  seen  in  Figure  4.3.  This  refined  mesh 
also  predicted  the  crack  growth  as  a  function  of  load  as 
shown  in  Figure  4.4.  While  no  experimental  data  is  available 
to  confirm  these  predictions,  it  is  interesting  to  note  that 
each  increment  of  crack  growth  advanced  over  several  nodes. 

Element  orientation  also  plays  a  role  in  fracture  load 
prediction.  Figure  4.5  shows  one  example  of  this  effect. 
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This  difference  in  element  orientation  accounts  for  the 
slight  upswing  near  the  middle  of  Figure  4.2. 

Finally,  the  program  shows  the  material  state  of  each 
element  as  the  load  increases.  This  data  is  plotted  at  se¬ 
lected  load  intervals  in  Figure  4.6.  Part  (a)  shows  the 
initial  formation  of  the  plastic  zone  at  a  load  of  10,000 
pounds.  As  the  load  increases  to  20,000  pounds  (Part  (b)), 
the  region  of  the  specimen  with  properties  in  the  first 
plasticity  section  increases  and  a  small  region  in  the 
second  and  third  sections  begin  to  form  at  the  crack  tip. 
Part  (c)  shows  the  expansion  of  all  three  regions  just  prior 
to  initial  fracture.  One  quadrant  of  the  specimen  is  also 
shown  in  (c)  to  indicate  the  relative  size  of  the  plasticity 
zones.  Figure  4.6d  shows  the  plasticity  zones  after  a  sig¬ 
nificant  amount  of  stable  cracking.  Note  that  after  initial 
fracture  and  unloading,  an  element  may  go  directly  from  an 
elastic  response  into  any  of  the  plasticity  sections  depend¬ 
ing  on  its  previous  progress  along  the  stress-strain  curve 
(strain  hardening).  Note  also  that  the  plasticity  zone  is 
still  increasing  in  size  as  the  crack  advances.  In  Figures 
4.6e  and  4.6f  the  plasticity  zones  move  partially  outside 
the  magnified  area  of  the  crack  tip  with  the  region  shown 
in  Figure  4.6f  being  the  plasticity  zone  at  fracture. 

PANEL2  was  also  used  to  analyze  several  panel  meshes 
to  determine  the  effects  of  using  true  stress-strain  and 
updating  specimen  geometry  during  each  load  increment.  No 
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significant  improvement  was  noted;  therefore,  only  the  engi¬ 
neering  stress-strain  programs,  PANEL1  and  FRACTURE,  were 
used  in  the  subsequent  analyses. 

IV. 3  Cracked  Panel  Summary 

The  analysis  of  a  centrally  cracked  panel  under  mono- 
tonically  increasing  load  using  the  modified  finite  element 
program,  PANEL,  demonstrated  that  unstable  fracture  predic¬ 
tion  using  the  FEM  is  highly  dependent,  on  the  mesh  size 
at  the  crack  tip.  Additionally,  the  prediction  of  stable 
fracture  requires  that  a  suitably  refined  mesh  be  extended 
along  the  entire  length  of  the  anticipated  crack  growth. 

For  monotonically  increasing  load,  accurate  predictions 
can  be  made  using  the  engineering  stress-strain  relation  and 
initial  specimen  geometry. 


CHAPTER  V 


TENSILE  SPECIMEN  ANALYSIS 

V.l  Introduction 

As  mentioned  in  Chapter  IV,  the  fracture  prediction 
programs  require  stress-strain  data  all  the  way  to  ultimate 
load.  This  data  is  not  generally  published  for  the  high 
strain  range.  Even  if  it  were,  the  scatter  in  properties 
might  introduce  error  into  the  analysis  since  published 
stress  capabilities  are  normally  statistical  minimums . 
Therefore,  the  entire  stress-strain  curve  was  determined 
experimentally  for  the  2024-T3  sheet  from  which  experimental 
specimens  were  fabricated.  This  data  was  converted  to  sec¬ 
tional  modulus  and  Poisson's  ratio  which  were  then  used  in 
the  program.  Models  of  the  tensile  test  specimens  were  also 
run  in  FRACTURE  to  evaluate  the  effectiveness  of  the  pro¬ 
gram. 
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V.2  Experimental  Study 

The  stress-strain  relation  for  the  0.125  inch  thick, 
2024-T3  aluminum  sheet  used  in  the  fractured  specimen  analy¬ 
sis  (Chapter  VI)  was  obtained  from  uniaxial  tensile  tests 
on  specimens  whose  dimensions  are  shown  in  Figure  5.1. 

These  specimens  were  loaded  to  fracture  on  a  Riehle  test 
machine  which  provides  calibrated  load  data.  Pin  to  pin  de¬ 
flections  for  the  specimens  were  obtained  from  a  spring 
loaded  potentiometer  attached  to  the  pins.  These  deflec¬ 
tions  were  recorded  as  a  function  of  the  applied  load.  The 
resultant  load-deflection  curves  are  plotted  in  Figure  5.2. 

The  loads  (L)  were  converted  to  engineering  stress  (<r) 
by  dividing  by  the  original  cross  sectional  area  (A0) : 

o’  =  L  /  A0  (5.1) 

Only  pin  to  pin  deflection  data  was  obtained  to  avoid 
damage  to  instrumentation  when  specimens  were  loaded  to  cat¬ 
astrophic  failure.  It  was  therefore  necessary  to  adjust 
the  pin  to  pin  deflection  (Dpp)  to  a  gage  deflection  (Dg) . 
This  was  accomplished  by  selecting  a  two  inch  gage  length 
on  the  neck  section  and  assuming  that  outside  this  region 
the  material  remained  elastic.  An  elastic  finite  element 
program,  based  on  published  Young's  modulus  and  Poisson's 
ratio,  was  then  used  to  determine  the  relative  elastic 


deflections  between  the  pin  and  a  point  on  the  gage  boundary 
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as  a  function  of  load.  Calling  this  deflection  Dr,  the  gage 
deflection  is  approximated  by 

Dg(L)  =  Dpp(L)  -  2.0  Dr(L)  L  (5.2) 

The  factor  two  in  Eq.  (5.2)  results  from  there  being  two  pin 
to  gage  boundary  regions.  While,  as  will  be  shown  shortly, 
some  plastic  region  exists  outside  of  this  gage  region,  the 
error  is  considered  small.  This  results  in  a  strain  (C)  of 
the  form 

C  =  Dg/2.0  (5.3) 

or,  in  view  of  (5.2), 

G  =  (Dpp/2.0)  -  Dr  L  (5.4) 

The  resulting  stress-strain  curve  is  shown  as  the  solid  line 
in  Figure  5.3.  Note  that  the  strain  at  fracture  was  init¬ 
ially  determined  from  the  0.824  inch  width  specimen.  The 
stress-strain  curve  was  then  approximated  by  the  six  linear 
sections  as  shown  in  Figure  5.3. 

V.3  Finite  Element  Analysis 

The  finite  element  models  from  the  mesh  generator 
program  were  run  in  the  FRACTURE  program  to  test  its  ability 
to  duplicate  the  load  deflection  curves  which  generated  the 
stress-strain  data  used  in  the  program.  Figure  5.4a  shows 
the  original  coarse  mesh  and  Figure  5.4b  shows  a  medium 
mesh.  The  medium  mesh  is  refined  in  the  area  of  the  fillet 
and  pin  sections.  The  fine  mesh  is  used  at  the  midsection 
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in  both  models  in  anticipation  of  fracture  in  this  area. 
Again,  only  one  quadrant  of  the  specimen  was  modeled  due 
to  symmetry. 

Figure  5.5  shows  a  comparison  of  the  results  of  the 
FEM  analysis  with  the  experimental  loads  and  deflections. 

For  the  0.824  inch  wide  specimen  (used  to  obtain  the  stress- 
strain  curve  for  the  program) ,  the  results  obtained  from 
both  meshes  are  very  accurate;  however,  the  medium  mesh, 
with  refinements  in  the  pin  and  fillet  areas,  gives  slightly 
more  accurate  results  in  the  elastic  and  fracture  regions. 

A  further  refinement  of  the  mesh  in  the  neck  area  was  found 
to  have  negligible  effect  on  the  results.  The  results  ob¬ 
tained  for  0.759  inch  wide  coarse  model  shows  good  agree¬ 
ment  with  experimental  results  except  at  fracture,  even 
though  the  material  properties  were  obtained  from  the  wider 
specimen.  When  a  medium  mesh  (results  not  shown  for  clarity) 
was  run  for  the  0.759"  model,  fracture  occurred  at  approxi¬ 
mately  the  sime  deflection  as  the  0.824"  model,  as  opposed 
to  the  larger  deflection  of  the  test  specimen.  This  sug¬ 
gests  that,  if  the  FEM  analysis  is  assumed  to  be  correct, 
the  difference  in  deflections  at  fracture  for  the  two  ex¬ 
perimental  specimens  is  not  accurate;  indeed,  the  difference 
was  traced  to  a  slight  anomaly  in  width  of  the  0.824"  speci¬ 
men.  To  compensate,  the  0.824"  load-deflection  curve  was 
extrapolated  out  to  the  0.759"  deflection  at  failure,  and 
an  adjustment,  shown  in  Figure  5.3,  was  made  to  the  program's 
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5.5.  FEM  Load  Deflection  Predictions  for  Tensile  Test 
Specimens . 
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stress-strain  data.  The  medium  mesh  was  then  used  in  the 
program  with  the  new  yield  region  5  material  properties. 
Correlation  with  the  experimental  results  was  much  better. 

The  remaining  difference  is  most  probably  due  to  the  fact 
that  the  deflection  at  fracture  should  be  higher  for  the 
0.824"  specimen  than  for  the  0.759"  specimen  since  there  is 
additional  deflection  in  the  elastic  region  for  the  0.824" 
specimen  due  to  its  increased  load  carrying  capability  at 
fracture.  The  adjusted  stress-strain  curve  was  then  used 
for  the  fracture  studies  of  Chapter  VI . 

Not  only  did  the  program  demonstrate  the  ability  to 
accurately  predict  the  load-deflection  curves  which  generated 
its  material  properties,  it  also  yielded  the  following  impor¬ 
tant  and  useful  data:  First,  note  on  Figure  5.5  that  local 
yielding  occurred  well  before  yielding  became  apparent  in 
the  load  deflection  curve.  Second,  the  program  provided 
data  on  the  progression  of  yield  through  the  specimen.  Fig¬ 
ure  5.6  shows  the  smoothed  yield  response.  Initial  yielding 
occurs  at  the  fillet  as  shown  in  Figure  5.6a.  At  a  load  of 
approximately  1000  pounds  below  the  apparent  yield,  this 
region  spreads  through  the  fillet  area  and  begins  at  the 
edge  of  the  pin  (while  the  yielding  at  the  pin  was  con¬ 
firmed  by  measuring  the  hole  after  fracture,  the  coarse¬ 
ness  of  the  mesh  in  this  region  may  not  have  given  an  accu¬ 
rate  map  of  the  yield  zone).  In  Figure  5.6c,  one  element 
at  the  fillet  has  moved  into  the  second  yield  region,  the 
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Figure  5.6  (cont).  Yield  Regions,  Tensile  Test  Specimen. 
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fillet  yield  region  has  expanded,  and  yielding  has  begun  on 
the  center  axis.  In  Figure  5.6d,  the  neck  area  is  almost 
entirely  involved  in  the  first  yield  region.  The  two  large 
white  areas  in  the  neck  are  still  elastic.  In  Figure  5.6e, 
the  entire  neck  is  involved  in  the  first  yield  region,  the 
second  yield  region  at  the  fillet  is  expanding,  and  second 
yield  has  occurred  at  the  center  axis.  The  second  region 
at  the  fillet  is  expanding  in  Fig.  5.6f,  g  and  h  with  the 
first  yield  region  moving  up  to  the  base  of  the  neck  while 
the  first  yield  zone  spreads  at  the  pin.  In  Figure  5.6i, 
third  yield  (dark  area)  is  progressing  in  much  the  same 
way  as  the  second  did  in  Figure  5.6f.  Also  in  Figure  5.6i, 
the  second  yield  region  is  entered  at  the  pin.  Figure  5.6j 
and  k  show  the  further  expansion  of  the  third  yield  zone  with 
Figure  5.6(1)  showing  the  beginning  of  the  fourth  yield  re¬ 
gion.  The  first,  second,  and  third  yield  zones  are  com¬ 
pressed  toward  the  base  of  the  neck  as  zone  four  expands  in 
Figure  5.6m  and  n,  with  the  pin  zone  continuing  to  expand. 
Fifth  yield  initiates  from  the  center  axis  as  shown  in 
Figure  5.6(o).  In  Figure  5 . 6p  the  pin  zone  increases  fur¬ 
ther  along  with  the  fifth  yield  region  while  zones  one,  two, 
three  and  four  are  pushed  further  toward  the  base  of  the 
neck.  Initial  fracture  occurs  in  Figure  5.6q  initiating 
from  the  center  axis  and  propagating  unstablely  to  the  edge 
as  shown  in  Figure  5.6r.  Note  that  the  FEM  prediction  for 
fracture  load  is  99.5%  of  the  experimentally  obtained  load. 
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Third,  as  mentioned  above,  the  FEM  analysis  predicted  an 
initiation  and  unstable  propagation  of  the  crack  from  a 
point  on  the  longitudinal  centerline  of  the  specimen.  This 
phenomenon  was  confirmed  on  the  test  specimens  by  placing 
the  fractured  surfaces  together  and  observing  that  the  end 
sections  fit  together  while  the  center  sections  did  not. 

This  is  due  to  the  increased  plastic  strain  on  the  outer 
sections  after  the  center  section  failed  and  unloaded.  Fail¬ 
ure  of  mildly  notched  tensile  specimens  from  the  center  axis 
has  also  been  reported  by  Drucker.^  Since  initial  yield 
occurred  at  the  fillet,  and  fracture  initiated  on  the  center 
axis,  it  can  readily  be  seen  that  fracture  initiation  loca¬ 
tion  can  not  be  predicted  by  using  the  maximum  elastic  stress 
location.  Fourth,  while  the  mesh  was  refined  on  the  mid¬ 
section  in  anticipation  of  failure  along  the  centerline, 
the  FEM  prediction  showed  that  fracture  occurred  off  the  mid¬ 
section  centerline  as  shown  in  Figure  5.6r.  Each  of  the 
tensile  test  specimens  also  broke  along  a  line  off  the  cen¬ 
terline.  Figure  5.7  shows  the  location  of  fracture  predicted 
by  the  FEM  analysis  and  as  occured  in  the  tensile  tests. 

The  test  specimens  failed  on  a  45°  line  through  the  thick¬ 
ness.  This  is  the  scatter  band  shown  in  the  figure  with 
the  experimental  location  shown  in  both  the  deflected  and 
undeflected  geometries.  Note  that  the  FEM  prediction  indi¬ 
cates  that  a  perfect  specimen  would  break  into  three  pieces. 
As  can  be  seen,  the  FEM  accurately  predicts  the  failure  loca¬ 


tion  . 


63 

Finally,  the  load  and  deflection  at  fracture  for  both 
the  experimental  test  and  the  FEM  analysis  are  shown  in 
Table  5.1.  The  load  predictions  for  the  medium  mesh  models 
are  extremely  accurate  (less  than  three  percent  for  the 
0.759"  specimen  and  less  than  one  percent  for  the  0.824" 
specimen) ,  and  deflections  are  also  a  good  approximation  of 
measured  values  (about  four  percent  and  one  percent  for  the 
two  respective  specimens).  As  previously  discussed,  accu¬ 
rate  prediction  of  fracture  deflection  is  highly  dependent 
on  accurate  material  maximum  strain  data. 

V.4  Tensile  Test  Analysis  Summary 

The  finite  element  program,  FRACTURE,  demonstrated 
the  following  capabilities  for  the  analysis  of  two  tensile 
test  specimens: 

1.  Ability  to  predict  load-deflection  curves, 

2.  Ability  to  demonstrate  the  importance  of  local 
material  properties, 

3.  Ability  to  provide  data  on  the  complete  field  re¬ 
sponse  for  the  specimen,  thus  a  better  understanding  of  the 
failure  process, 

4.  Ability  to  predict  fracture  initiation  location, 
both  with  respect  to  the  midsection  and  longitudinal  axis, 
and 

5.  Ability  to  predict  load  and  deflection  at  frac¬ 


ture  . 


Locations,  FEM  Predictions  and  Experimental  Results  for  Tensile 


TABLE  5.1 

COMPARISON  OF  EXPERIMENTAL  AND  FEM  LOADS 
AND  DEFLECTION  AT  FRACTURE 


SPECIMEN 


0.759" 


0.824" 


LOAD 

(lb) 


DEFL. 

(in) 


LOAD 

(lb) 


DEFL. 


EXPERIMENTAL 


7180 


0.3750 


7580 


0.3125 


COARSE 

MESH 


6953 


0.2928 


7532 


0.3035 


PERCENT 

ERROR 


3.16 


21.92 


0.63 


2.88 


MEDIUM 

MESH 


7015 


0.3588 


7539 


0.3158 


PERCENT 

ERROR 


2.30 


4.32 


0.54 


1.06 
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CHAPTER  VI 

CRACKED  SPECIMEN  ANALYSIS 

VI .  1  Introduction 

The  tensile  test  specimen  analysis  demonstrated  the 
effectiveness  of  the  current  method  in  predicting  experimen¬ 
tal  behavior  for  mild  stress  concentrations  and  the  panel 
study  showed  that  for  sufficiently  fine  meshes,  this  method 
can  also  predict  load  at  fracture  in  specimens  with  severe 
stress  concentrations  (cracks).  Unfortunately,  no  experi¬ 
mental  deflection  data  was  presented  for  the  panel  study; 
therefore,  the  ability  of  the  program  to  predict  load  de¬ 
flection  curves  for  severely  notched  specimens  could  not  be 
addressed  without  further  testing.  To  obtain  the  needed 
data,  three  tensile  test  specimens  were  notched  and  loaded 
to  fracture.  Comparisons  of  this  experimental  data  and  the 
finite  element  predictions  were  made  and  are  presented  in 
this  chapter. 
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VI .2  Experimental  Tests 

The  testing  procedure  for  notched  specimens  was  the 
same  as  that  presented  in  Chapter  V,  and  the  three  specimens 
tested  were  made  from  blanks  of  the  same  dimensions  as  the 
0.824"  wide  specimen  of  Chapter  V.  Sharp  notches  of  lengths 
0.008",  0.023",  and  0.129"  were  then  introduced  in  the  blanks 
on  one  edge  of  the  neck  at  the  centerline.  The  0.003"  and 
0.023"  cracks  were  obtained  using  an  X-Acto  knife  biode, 
and  the  0.129"  crack  was  machined  on  a  band  saw  with  the 
final  tip  also  being  formed  by  an  X-Acto  knife. 

VI .3  Finite  Element  Analysis 

Introduction  of  the  single  edge  notch  (SEN)  in  the 
specimen  removed  one  plane  of  symmetry  necessitating  the 
use  of  a  two  quadrant  finite  element  model.  The  overall 
mesh  is  shown  in  Figure  6.1a  with  details  for  the  different 
crack  lengths  shown  in  Figure  6.1b,  c,  and  d. 

Figure  6.2  shows  the  same  type  of  result  that  was  dem¬ 
onstrated  in  Chapter  IV,  Figure  4.2;  that  is,  the  accuracy 
of  the  FEM  predictions  of  fracture  load  is  highly  dependent 
on  the  element  size  at  the  tip  of  the  crack. 

The  load  deflection  curves  for  two  different  size 
elements  are  shown  in  Figure  6.3.  Note  that  the  shape  of 
the  predicted  curves  are  essentially  the  same.  With  large 
elements  at  the  tip,  the  CST  elements  can  not  model  the  large 
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gradients  at  the  tip,  and  therefore,  the  load  and  deflec¬ 
tion  exceed  the  experimental  values.  The  slight  upward 
trend  in  the  load  for  the  finer  mesh  is  due  to  the  size  of 
the  minimum  load  increments  used  to  reduce  the  computational 
time.  As  the  mesh  is  refined,  the  computational  time  in¬ 
creases  due  to  increased  band  width  of  the  stiffness  matrix 
and  the  increased  number  of  nodes  and  elements.  The  compu¬ 
tational  time  can  be  reduced  by  increasing  the  minimum  load 
increment;  however,  the  predicted  load  and  deflection  at 
fracture  are  affected  since  elements  tend  to  remain  stiffer 
during  the  loading  process. 

The  load  deflection  curves  for  the  three  different 
size  notches  are  shown  in  Figure  6.4.  The  results  improve 
as  the  crack  size  increases.  This  can  be  attributed  to  in¬ 
creasing  the  minimum  load  increment  to  allow  enough  time  to 
advance  the  crack  along  the  additional  specimen  width  for 
smaller  cracks.  Also  the  deflections  are  slightly  low  for 
each  given  load.  This  same  effect  can  be  observed  in  Fig¬ 
ure  5.5  for  the  tensile  specimen  coarse  mesh.  As  discussed 
in  Chapter  V,  the  tensile  predictions  were  improved  by  re¬ 
fining  the  mesh  in  the  area  of  the  fillet  and  pin. 

Finally,  the  yield  regions  for  a  0.023"  initial  crack 
size  are  plotted  for  selected  loads  in  Fibure  6.5.  The 
initial  yield  zone  formation  is  shown  in  Figure  6,5a  and  b. 
Figure  6.5c  and  d  show  the  first,  second,  third  and  fourth 
yield  regions  expanding  outward  from  the  crack  tip.  At  the 
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Figure  6.1.  JEM  Models  of  Cracked  Specimens. 


(d)  Detail  of  Mesh  for  0.129  In.  Crack 


Figure  6.1  (cont).  FEM  Models  of  Cracked  Specimens 
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.3.  Effects  of  Load  Increment  and  Element  Size 
for  the  Cracked  Specimen  Analysis. 


Figure  6.4.  Load  Deflection  Curves  for  Cracked  Specimen  Analysis. 
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load  that  corresponds  to  Fig.  6.5c,  yielding  also  begins  at 
the  hole.  The  first  yield  zone  extends  past  the  enlarged 
area  of  the  figure  (as  shown  in  Figure  6.5e)  with  the  entire 
specimen  shown  in  Figure  6.5f  for  the  same  load.  Figures 
6.5g,  h  and  i  show  the  further  advance  of  the  yield  zones 
as  load  increases.  At  load  level  corresponding  to  Figure 
6.5g,  all  five  zones  are  present.  In  Figure  6.5h,  the  speci¬ 
men  has  reached  the  load  at  which  initial  fracture  occurs  and 
Figure  6.5i  shows  the  zones  just  prior  to  unstable  fracture. 
Notice  that  the  entire  neck,  except  at  the  tip,  is  in  the 
second  yield  region. 

VI .4  Cracked  Specimen  Analysis  Summary 

The  FEM  analysis  of  the  sharply  notched  specimens  des¬ 
cribed  in  this  chapter  further  verifies  the  abilities  demon¬ 
strated  in  the  previous  two  chapters,  and  confirms  the  abil¬ 
ity  to  predict  specimen  deflections  for  sharply  notched 
specimens . 


CHAPTER  VII 


CONCLUSIONS  AND  RECOMMENDATIONS 


The  finite  element  program  developed  during  this  study 
has  been  shown  to  predict  load  deflection  curves,  load  and 
deflection  at  fracture,  fracture  paths,  initiation  sites, 
crack  growth,  and  stable  or  unstable  crack  propagation.  The 
accuracy  of  the  method  is  highly  dependent  on  the  element 
size  along  the  crack  path.  The  use  of  properly  refined 
meshes  yields  accurate  results. 

The  method  is  completely  general  in  that  it  can  anal¬ 
yze  any  two-dimensional  isotropic  structure  subjected  to 
plane  stress  and  uniaxial  loading.  The  loading  restriction 
can  be  removed  by  substituting  a  failure  criterion  which  is 
more  suitable  than  the  maximum  strain  criterion.  While  this 
method  represents  a  valuable  design  tool  for  a  limited  class 
of  problems,  it  more  importantly  demonstrates  the  potential 
of  the  finite  element  method  for  the  direct  prediction  of 
fracture.  The  current  program  required  only  minor  altera¬ 
tions  to  a  standard  plane  stress  finite  element  program  to 
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give  accurate  analysis  of  elasto-plastic  fracture  problems. 
Since  finite  element  programs  have  already  been  written  to 
analyze  plane  strain,  three  dimensional  and  thermal  loading 
problems,  if  similar  modifications  could  be  made  to  these 
programs,  then  it  would  be  possible  to  directly  predict 
fracture  for  these  cases  under  monotonically  increasing 
load.  The  only  restriction  on  this  approach  appears  to  be 
the  computer  storage  and  computational  time.  These  become 
less  significant  as  the  program  effectiveness  improves,  re¬ 
finements  such  as  substructuring  are  incorporated,  better 
solution  techniques  are  found,  and  as  computer  capabilities 
continue  to  expand. 

The  direct  prediction  of  cyclic  fracture  would  be  an 
even  more  valuable  application  of  the  approach  contained  in 
this  work.  Again  the  basic  procedures  developed  in  this 
study  should  apply  with  appropriate  modifications. 

Each  of  these  capabilities  needs  to  be  verified,  but 
the  excellent  results  obtained  in  the  present  study  suggest 
that  the  concept  is  valid  and  worthy  of  further  development. 
The  rewards  for  such  a  work  could  be  enormous . 
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APPENDIX  I 


THE  FINITE  ELEMENT  PROGRAM,  FRACTURE 


This  appendix  contains  the  finite  element  FRACTURE. 

A  sample  of  the  program  output  is  also  included  at  the  end 
of  the  program. 
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